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Development  of  a  new  boundary  layer  code  has  reached  the  status  to  permit  meaningful 
applications.  Computations  have  been  carried  out  for  the  initiation  of  separation  in  the  symmetry 
plane  of  a  prolate  spheroid  of  slenderness  ratio  1/4,  impulsively  started  into  forward  motion  at  zero 
incidence  and  also  at  50  degrees  angle.  This  case  serves  as  validation  by  comparing  with 
previously  published  results  of  Xu  and  Wang  (ref.  1),  and  also  demonstrates  that  our  method  gives 
informaron  of  flow  near  the  rear  stagnation  point  not  available  in  the  literature.  More  studies  have 
been  performed  on  the  optimization  of  surface  suction  to  delay  or  prevent  the  unsteady  separation 
for  an  impulsively  started  circular  cylinder.  Here  the  methodology  should  be  of  interest.  Work  on 
an  unsteady  three-dimensional  thin-layer  Navier-Stokes  code,  however,  is  progressing  slowly. 
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Technical  Rept  for  AFOSR  88-0229,  for  the  period  Oct.  '89  to  Feb.  '90 
by  S.F.  Shen,  Cornell  University 

We  have  reached  the  termination  date  of  the  Grant  88-  0229,  except  for  a  new 
minimal  extension  to  enable  us  to  bring  our  major  effort,  Wu's  study  of  the  unsteady 
separation  process  on  a  three-dimensional  body,  to  at  least  some  concrete  conclusion  by  the 
end  of  August,  1990.  While  we  are  grateful  for  Capt.  Helin's  sustained  interest  in  our 
project,  the  amount  available  to  us  from  AFOSR  will  not  cover  half  the  cost.  We  shall  try 
to  manage  with  forced  reductions  of  expenditure  in  the  last  few  months  and  the  principal 
investigator  is  taking  zero  compensation. 

The  progress  of  the  project  between  Oct.  '89  and  Feb.  '90  in  all  three  active  studies 
outlined  in  the  last  Technical  Report  ( Oct.  ’89)  is  summarized  as  follows.  They  consist  of: 

1)  Three-dimensional  unsteady  separation  over  an  impulsively  started  prolate 
spheroid,  at  zero  and  also  50  degrees  of  angle  of  attack; 

2)  Feasibility  study  of  controlling  two-dimensional  unsteady  separation  over  a 
circular  cylinder  by  surface  suction; 

3)  Development  of  an  unsteady  thin-layer  ('parabolized')  Navier-Stokes  code  for 
boundary-layer  behavior  subsequent  to  the  initiation  of  separation. 

Project  1)  has  been  carried  out  by  T.Wu.  During  this  period  the  focus  has  been  a 
NEW  three-dimensional  unsteady  boundary-layer  Eulerian  code  for  flow  in  the  symmetrical 
plane.  After  careful  validation  it  is  now  in  operational  form.  The  first  application  was  to 
recalculate  the  case  in  the  published  work  of  Xu  and  Wang  (ref.  1),  having  to  do  with  the 
process  leading  to  initial  separation  in  the  plane  of  symmetry.  Our  procedure  in  principle 
has  two  major  advantages:  the  multizone  formulation  eliminates  the  difficulty  associated 
with  the  geometrical  singularity,  and  the  initial-value  approach  allows  accurate  treatment  of 
the  moving  stagnation  points.  A  paper  has  been  prepared  for  presentation  at  the 
International  Symposium  on  Nonsteady  Flows  at  the  joint  ASME  and  CSME  (  Canadian) 
meeting,  June  '90  in  Toronto,  Ont.  A  copy  is  attached  as  Appendix  A. 

Project  2)  has  been  the  assignment  of  Z.  Xiao.  The  basic  idea  is  to  exploit  our 
time-accurate  Eulerian  unsteady  boundary-layer  code,  to  study  the  effects  on  the  initiation 
of  unsteady  separation  due  to  various  parameters  related  to  the  variable  motion,  as  well  as 
the  boundary  conditions  that  may  suggest  strategies  for  separation  control.  The  Eulerian 
code  must  be  able  to  locate  incipient  separation  that  has  been  made  clear  by  our  persistent 
demonstration  via  the  Lagrangian  interpretation.  In  the  Eulerian  description,  most  experts 
now  generally  agree  that  the  precursor,  or  signature,  of  imminent  separation  is  a  locally 
very  sharp  'spike'  in  the  displacement  thickness.  We  first  establish  that  in  the  bench-mark 
case  of  an  impulsively  started  circular  cylinder,  the  van  Dommelen  and  Shen  (ref.2) 
solution,  recast  as  the  displacement  thickness  along  the  wall,  is  indeed  duplicated  by  our 
new  Eulerian  code.  Preliminary  results  have  been  obtained  on  the  efficacy  of  delaying  or 
eliminating  the  occurrence  of  the  'spike',  hence  separation,  with  the  applicadon  of  various 
surface  suction  arrangements.  These  were  included  in  an  invited  presentation,  by 
S.F.Shen,  at  the  Workshop  on  Analytical  Methods  in  Unsteady  Separation,  sponsored  by 
Army  ARO,  January  '90,  Columbus  OH.  A  copy  of  the  handout  of  the  talk  was  submitted 
with  our  Technical  Report  for  Oct.  '89. 


We  hope  for  more  results  on  optimization  of  controlling  unsteady  separation  before 
his  departure.  Because  of  the  budget  limitation,  however,  Mr.  Xiao's  association  with 
AFOSR-sponsored  research  will  be  terminated  no  later  than  March  1,  1990. 

Project  3)  has  occupied  the  attention  of  Dr.  J.  S.Kim  since  his  arrival  from  S. 
Korea  in  may  '89.  As  described  in  earlier  Technical  Reports,  Dr.  Kim's  attack  follows  his 
previous  experiences  with  the  Navier-Stokes  solver.  To  rewrite  some  of  the  program  for 
use  in  the  Cornell  Supercomputer  facility  proved  to  be  quite  time  consuming,  esp.  for 
applications  targeted  for  three-dimensional  bodies.  Considerable  debugging  efforts  have 
been  necessary  during  the  present  period,  and  so  far  the  test  trials  involve  only  simple 
geometry  and  small  numbers  of  time  steps.  As  Dr.  Kim  is  obligated  to  return  to  Korea  on 
April  1,  1990,  and  no  funding  for  possible  replacement,  it  seems  unrealistic  to  expect  the 
final  achievement  of  this  work  to  be  more  than  preliminary  exploration.  A  description  of 
Kim’s  work  can  be  found  in  Appendix  B. 

With  regard  to  publication  and  presentation,  the  paper  by  Shen  and  Wu  (AIAA 
Paper  88-  3542  )  is  now  'in  the  press'  of  the  AIAA  Journal,  after  a  prolonged  unfortunate 
status  as  a  'missing  manuscript'  in  the  editor's  office.  In  January  S.F.Shen  gave  a 
colloquium  on  'Unsteady  Boundary  Separation  from  the  Lagrangian  View'  at  the 
University  of  Arizona,  Tuscon,  AZ,  followed  by  the  invited  talk  of  the  same  title  (but 
somewhat  different  content,  including  some  of  the  separation  control  results)  at  the 
Columbus  Workshop.  Wu's  results  to-date  will  be  presented  at  the  ASME  meeting  in 
Toronto  this  June,  already  mentioned  above. 
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B .  Code  Development  for  the  Unsteady  Thin  Layer  Navier-Stokes. 
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A  MULTI-ZONE  TIME-MARCHING  TECHNIQUE  FOR  UNSTEADY 
SEPARATING  THREE-DIMENSIONAL  BOUNDARY-LAYERS  AND  ITS 
APPLICATION  TO  THE  SYMMETRY-PLANE  SOLUTION  OF  AN 
IMPULSIVELY-STARTED  PROLATE  SPHEROID 

Tzuytn  Wu*  and  Shan-Fu  Shen** 

Sibley  School  of  Mechanical  and  Aerospace  Engineering 
Cornell  University 
Ithaca,  New  York,  14853 


ABSTRACT 

Recent  interest  in  unsteady  separation  and  separated  flows  brings 
up  the  need  of  an  accurate  and  efficient  computational  scheme  for 
general  unsteady  three-dimensional  boundary-layer  flows.  Resolution 
of  the  singular  behavior  at  separation  is  a  delicate  problem.  The  task 
is  further  complicated  by  the  geometrical  singularity  and  the  non¬ 
stationary  stagnation  point.  The  present  paper  proposes  a  numerical 
scheme  to  sidestep  these  difficulties.  As  the  first  stage  of 
development,  the  simpler  problem  of  the  symmetry-plane  solution  of 
the  laminar  boundary-layer  over  an  impulsively  started  prolate 
spheroid  is  calculated.  Results  show  that  the  present  Eulerian 
calculation  satisfactorily  captures  the  singular  behavior  of  the 
boundary-layer  when  separation  is  approached.  Comparison  with  Xu 
and  Wang's  recent  results  and  those  for  the  two-dimensional  elliptic 
cylinder  calculated  by  the  Lagrangian  method  are  also  made. 
Discussions  of  the  results  for  unsteady  separation  at  zero,  small  and 
large  incidences  are  presented. 
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difference  operator  in  finite-difference  equations 
indicating  difference:  also,  displacement  thickness 
stretching  factor  of  the  mapping  in  a-direction 
average  operator  in  finite-difference  equations 
density 
wall  shear 


Subscripts 
o  stagnation 

s  separation 

x,y,z  components  in  x,  y  and  z  directions  respectively 

,  indicating  derivatives 

Superscripts 

n  n1*1  time  level  in  finite-difference  equations 

O  quantities  in  nose  region 

O  quantities  in  tail  region 
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mapping  constants  in  [3-direction 
eccentricity  of  spheroid 
metric  coefficients 
pressure 

transformed  coordinate  defined  in  equation  (4) 
time 

thickness  ratio 

velocity  component  in  x-direction 
value  of  u  at  edge  of  boundary-layer 
velocity  component  in  y-direction 
value  of  v  at  edge  of  boundary-layer 
derivative  of  v  in  y-direction 
value  of  vv  at  edge  of  boundary-layer 
constants  in  the  velocity  distribution  at  edge  of 
boundary -layer 

velocity  component  in  z-direction 

streamwise  coordinate  along  the  major  axis  of  prolate 

spheroid 

circumferential  coordinate 

boundary-layer  coordinate  normal  to  x  and  y  direction 
angle  of  attack 
computational  coordinates 
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The  separation  of  three-dimensional  laminar  boundary-layer  has 
long  been  an  intriguing  subject  of  both  basic  and  practical  interest. 
Even  in  the  steady  case,  the  separation  "pattern’'  on  the  body  surface, 
as  observed  in  the  experiments,  takes  a  variety  of  forms.  Recent 
progress  in  the  qualitative  theory  and  the  actual  computations  have 
greatly  improved  our  understanding  of  the  complex  phenomena  ([  1  - 
11]).  In  the  steady  flow,  zero  wall  shear  signifies  separation  and  the 
breakdown  of  the  boundary-layer  equation.  The  usual  analysis  of 
steady  separation  is  to  focus  on  the  limiting  streamlines  at  the  wall. 
For  the  unsteady  case,  it  is  now  well-known  that  separation  starts 
"off'  the  wall  and  neither  the  wall  shear  nor  the  wall  streamlines  can 
serve  as  reliable  indicators  of  separation.  But  very  little  is  known 
about  the  details  of  the  phenomenon.  Even  numerical  examples  are 
hard  to  come  by,  since  efficient  and  reliable  computational  schemes 
for  calculating  the  general  unsteady  three-dimensional  boundary  layer 
flow  still  are  not  well  developed. 

In  steady  flow,  the  separation  of  the  boundary-layer  over  a  prolate 
spheroid  at  various  incidence  angles  has  been  the  subject  of  several 
computational  studies  due  to  Wang  ([6-8],  [12-14]).  Two  different 
types  of  separation,  "open"  and  "closed",  based  on  the  wall  limiting 
streamlines  analysis  are  proposed.  In  his  terminology,  the 
convergence  of  the  limiting  streamlines  defines  the  separation  line  on 
the  body.  At  small  incidence,  a  "closed"  type  separation  prevails 
where  the  separation  line  completely  separates  the  limiting  streamlines 
originating  from  different  stagnation  points.  At  moderate  incidence, 
the  separation  line  is  approached  on  two  sides  by  limiting  site  anilines 
originating  from  the  same  stagnation  point.  This  tvpr  of  separation  is 
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then  classified  as  "open''  separation,  since  two  separation  lines,  one 
on  each  side  of  the  spheroid  are  formed  and  are  not  connected  on  the 
lee-side  symmetry  plane.  In  a  recent  paper  of  Xu  and  Wang  [151,  the 
evolution  of  the  skin  friction  on  the  symmetry  plane  of  an  impulsively 
started  spheroid  is  studied.  They  suggest  an  unsteady  analogy  of  the 
above  "open"  and  "closed"  separations,  the  term  "separation"  used  in 
their  study  referring  to  vanishing  skin  friction  in  the  streamwise 
direction.  An  off-svmmetry-plane  calculation  of  the  above  problem 
has  been  performed  by  Ragab  (161.  Focus  of  his  work  is  on  the 
crosswise  separation  at  high  incidence.  The  computational  domain 
covers  only  the  front  quarter  of  the  prolate  spheroid.  His  results  of 
skin  friction  lines  support  the  "open"  type  separation  advocated  by 
Wang  for  unsteady  three-dimensional  (lows.  Following  van 
Dommelen  and  Shen  (17,  18],  unsteady  boundary-layer  separation  is 
better  defined  by  the  termination  of  the  boundary-layer  solution  in  a 
singularity.  In  the  two-dimensional  problem,  the  occurrence  of  a 
singularity  and  its  structure  have  been  well  established  (e.g.  (17-201). 
In  the  three-dimensional  case,  this  singularity  may  form  a  curve 
whose  projection  on  the  wall  can  then  be  properly  defined  as  the 
"separation  line"  (Shen  (21)).  The  concept  is  consistent  with  that  of 
van  Dommelen  and  Cowley  [22|,  in  which  the  authors  examine  the 
asymptotic  structure  of  the  unsteady  three-dimensional  separation  in  a 
Lagrangian  coordinates  system.  In  fact,  the  structure  of  separation  is 
found  to  be  quasi  two-dimensional. 

Confirmation  of  this  singularity  at  separation  was  not  part  of  Xu 
and  Wang's  calculation,  but  is  the  focal  point  of  the  present  paper.  In 
order  to  do  so,  the  numerical  scheme  has  been  modified.  It  is  noticed, 
for  instance,  that  no  details  of  the  flow  near  the  trailing  edee  were 
included  in  their  results.  As  pointed  out  in  Cebeci  et  al.  ([10|,  [1 1], 
(23|),  a  "geometrical  singularity"  is  present  at  both  ends  of  the  prolate 
spheroid  when  a  body-oriented  coordinates  system  is  used.  In  Xu 
and  Wang  (15),  no  provision  was  made  for  the  singularities,  and  their 
method  might  lead  to  local  difficulty  and  inaccuracy,  particularly  for  a 
slender  spheroid.  Another  problem  of  the  space-marching  schemes 
usually  adopted  in  solving  a  boundary-layer  type  equation  is 
associated  with  the  stability  consideration.  In  the  unsteady  case,  the 
grid  size  in  the  reversed  flow  region  is  limited  by  the  domain-of- 
dependence  rule  (CFL  criterion).  In  particular,  the  restriction  on  the 
time  step  can  become  so  stringent  that  the  calculation  is  practically 
prevented  from  proceeding  further  [24],  Still  another  technical 
difficulty  of  the  space-marching  method  is  the  requirement  of  a 
starting  profile  to  initiate  the  spatial  integration.  This  is  not  available 
in  the  general  case  involving  arbitrary  body  motion  for  which  the 
stagnation  point  is  non-stadonary. 

The  present  paper  is  actually  the  first  phase  of  a  more  ambitious 
research  program  dealing  with  general  three-dimensional  unsteady 
separation  over  an  arbitrarily  moving  body.  In  developing  a  solution 
scheme,  we  abandon  the  conventional  space-marching  technique  and 
choose  instead  a  more  classical  initial-value  approach.  That  is,  the 
solution  is  advanced  in  the  time  direction  only.  This  allows  the 
calculation  of  more  general  unsteady  flows  as  it  does  not  require  a 
fixed  stagnation  point  to  construct  the  flow  downstream.  In 
conjunction  with  an  implicit  time-marching  scheme,  a  flow  dependent, 
one-sided  spatial  difference  is  applied  to  discretize  the  convection 
terms  in  accordance  with  the  CFL  criterion.  The  severe  restriction  on 
the  time  step  is  thus  relieved.  Finally,  to  overcome  the  difficulty  of 
the  aforementioned  "geometrical  singularity",  we  propose  a  multi- 
domain  attack.  Suitable  transformation  which  removes  the  singularity 
is  applied  at  both  front  and  tail  regions.  Calculations  are  performed 
separately  in  each  region,  and  the  solutions  are  patched  on  the 
interfaces  by  interpolation  procedure. 

For  an  orderly  development,  we  limit  ourselves  first  to  redo  the 
case  of  a  prolate  spheroid,  following  Xu  and  Wang  [15],  The 
following  presents  mainly  the  boundary- layer  development  in  the 
symmetrical  plane  of  an  impulsively-started  prolate  spheroid.  The 
governing  conations  and  the  relevant  coordinates  transformation  are 
given  in  the  next  section.  Details  of  the  numerical  scheme  and  the 
solution  procedure  are  also  described.  Calculations  are  performed  for 
three  different  angles  of  attack  :  0°,  6°  and  50°.  Comparison  with  Xu 
and  Wang  [15]  are  made.  In  addition,  the  results  are  compared  with 
those  for  the  two-dimensional  elliptic  cylinder  calculated  by  the 
Lagrangian  method.  The  significance  of  the  findings  of  the  first 
occurrence  of  three-dimensional  separation  on  the  lee-side  and  wind- 
side,  at  small  and  large  angles  of  incidence  are  discussed.  If  and 
when  and  how  different  unsteady  three-dimensional  separation 


patterns  are  generated  and  developed  can  only  be  answered  bv 
computation  of  the  flow  off  the  symmetry  plane.  Much  more  remains 
to  be  done. 

FORMULATION  OF  THE  PROBLEM 

The  geometry  and  the  boundary-layer  coordinate  system  of  a 
prolate  spheroid,  following  Xu  and  Wang  [15],  are  specified  in  Fig.l. 
The  metric  coefficients  hj,  hy  and  hj  are  defined  as 


hj  =  xj  *~e~  .  hs  =  V(l-e-)U-x') ,  li3  =  1  ;  e  -  l-a'ci- 

where  t/c  denotes  the  ratio  of  minor  to  major  axis  (axes  ratio).  After  a 
proper  scaling  with  the  free  stream  velocity  and  the  half-chord  length 
of  the  prolate  spheroid,  the  non-dimensional  forms  of  the  boundary- 
layer  equations  on  the  symmetry'  plane  are,  in  conventional  notation.” 
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Crosswise  momentum  equation  is  identically  zero  on  the  symmetry 
plane.  A  supplementary  equation  can  be  obtained  by  differentiating  it 
with  respect  to  the  circumferential  coordinate  y  : 


dvv  u  9vv  3vv  v2  uvv  9h->  -1  d2p  d2vv 

-^  +  RT^r  +  w”5F  +  hT  +  h^dr  =  ^:^  +  ^r 


(3) 


where  vv  stands  for  dv/dy.  Based  on  the  velocity  profiles  u  and  vv, 
the  skin  friction  and  the  displacement  thickness  are  defined  by 

PQ 
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o 

The  above  governing  equations  are  then  subject  to  a  Rayleigh  type 
initial  condition  for  an  impulsively-started  motion.  As  usual,  no  slip 
boundary  conditions  are  imposed  on  the  wall,  and  the  velocities  at  the 
outer  edge  of  the  boundary-layer  must  approach  the  values  specified 
by  the  inviscid  potential  flow  : 


u  =  vy  =  0  atz  =  0 

u  =  Ue  =  —  ~T2~U2 1  voO-x2)1/2  cos  &  +  v90  i  X  sin  a  cos  v  ] 


vy  =  Vyc  =  V9q  sin  a  cos  y 


where  &  =  angle  of  attack. 
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As  already  mentioned  in  the  introduction,  a  difficulty  of  the 
present  spherical  coordinates  system  is  that  the  metric  coefficients  hi, 
h2  become  singular  at  both  x=-l  and  x=+l.  A  transformation  which 
removes  this  geometric  singularity  was  proposed  by  Cebeci  et  a). 
[10],  [23],  Define 


2S  _  hydx  _  v/ 1  -e2x2~dx 
S  h2  Vu?(l-x2) 
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where  A  and  B  are  mapping  constants.  With  transformation  (41,  and 
further  let 


x  =  S  cos  y 
y  =  S  sin  y 
l  =  z 


u  =  u  cos  y  +  v  sin  y 
v  =  -u  sin  y  +  v  cos  y 
w  =  w 


Equations  1 1),  (2)  and  (3)  then  become 
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Equations  (5),  (6),  and  (7)  are  now  regular  at  x=-l.  Similar 
treatment  can  be  applied  to  remove  the  singularity  at  the  rear  end  x=+i 
and  we  shall  not  restate  the  derivation  here.  In  the  following  context, 
quantities  in  the  transformed  coordinates  near  ihe  trailing  edge  of  the 

prolate  spheroid  will  be  denoted  by  (  ). 

SOLUTION  METHOD 

The  solution  domain  is  decomposed  into  four  regions  as  shown  in 
Fig. 2.  The  transformed  coordinate  system  described  in  the  last 
section  is  used  in  the  nose  and  tail  regions  to  remove  the  geometry 
singularity.  Each  region  is  overlapped  with  the  adjacent  ones. 
Calculations  are  performed  separately  in  different  regions,  the  required 
boundary  conditions  on  the  interface  of  one  region  are  then  obtained 
from  interpolating  the  values  at  the  interior  points  of  the  other  region. 
Successful  results  where  similar  composite-grid  techniques  have  been 
applied  were  reported  in  several  previous  investigations  of  Navier- 
Stokes  flows  ([25],  [26]),  transonic  flows  ([27],  [28],  [29]),  Euler 
equations  ([29],  [30]),  and  the  combined  Navier-Stokes,  Euler  and 
boundary-layer  equations  [31].  The  principal  advantages  of 
composite-grid  are  the  generation  of  grids  for  regions  of  complicated 
geometry,  and  an  easier  implementation  of  different  formulations  in 
different  regions.  However,  the  "patching"  of  the  solutions  between 
different  regions  needs  attention.  Care  must  be  taken  so  that  the 
interpolation  procedure  on  the  interface  will  not  destroy  the  stability 
and  accuracy  of  the  interior  solution.  See,  e.g.  Chesshire  and 
Henshaw  [25].  In  the  following  calculations,  a  cubic  Lagrangian 
interpolation  formula  was  used.  No  apparent  error  and  instability 
were  observed  in  this  particular  numerical  example. 

The  physical  domain  is  more  conveniently  transformed  into  a  finite 

computational  region  (a, (5)  by  using  appropriate  mapping  functions. 
In  the  main  region,  for  example, 

x  =  (x2-X|)  [(l-e)(3-2a)  a2  +  ea)]  +  x,  0Sa<l 


0  <  (3  <  I 


in  the  a-direction  gives  the  mesh  clustering  effect  near  the  leading  and 
trailing  edges  where  the  pressure  gradient  varies  rapidly.  The  time¬ 
like  P  coordinate  removes  the  singular  behavior  in  the  flow  at 

impulsive  ttan-up.  A  uniform  gnd  is  drawn  in  the  a-(3  plane,  and  the 
governing  equations  are  discretized  in  terms  of  these  computational 
coordinates. 

In  contrast  to  the  conventional  space-marching  technique  in 
solving  the  boundary-layer  type  equation,  we  next  formulate  strictly 
an  initial-value  problem  in  which  the  flow  variables  in  each  region  are 
advanced  in  the  time-direction  only.  A  similar  approach  was  already 
carried  out  in  Van  Dalsem  and  Steger  [32]  who  used  a  time-relaxation 
scheme  in  then-  calculations.  We  seek  to  reduce  the  computation  effort 
at  each  time  step  by  adopting  a  non-iterative  predictor-corrector 
formulation,  Douglas  (33|.  TTiat  is,  in  advancing  the  time  step  from  r. 
to  n+1,  a  predicted  value  is  first  evaluated  at  n+1/2  where  all  the  non¬ 
linear  coefficients  in  the  momentum  equations  are  linearized  by  using 
the  values  frozen  from  the  previous  time  n.  Thus  the  governing 
equations  are  decoupled  in  the  predictor  step.  A  backward  Euler 
differencing  in  time  is  preferred  in  solving  the  predictor  value  for  the 
purpose  of  damping  out  the  numerical  error.  Second  order  accuracy 
of  temporal  difference  can  be  restored  by  following  a  Crank-Nicolson 
like  correction  step. 

In  both  the  predictor  and  the  corrector  steps,  a  flow  dependent, 
one-sided  spatial  difference  is  applied  to  the  convection  term  in 
accordance  with  its  hyperbolic  character  ibackward  difference  if  u>0 
and  forward  difference  for  u<0).  Unlike  the  usual  marching  scheme 
where  flow  quantities  needed  to  account  for  the  backflow  influence  are 
grabbed  from  the  previous  time  level.  By  appealing  to  the  CFL 
criterion,  this  explicit  formulation  will  then  restrain  the  time  step  used 
in  the  discretization  procedure.  This  is  particularly  important  as  the 
restriction  on  the  step  size  can  be  very  stringent  when  separation  is 
approached  [24].  In  the  present  scheme,  an  implicit  time-advancing 
approach  is  used,  and  upwind  differences  are  performed  at  the  current 
time  level  throughout  the  whole  domain.  The  CFL  condition  is 
satisfied  automatically,  allowing  us  to  relax  the  time-step  limitation. 

The  formulation  is  very  similar  in  both  the  predictor  and  corrector 
steps.  By  applying  the  above  discretization  procedures  to  (2)  and  (3). 
the  following  are  the  resulting  finite  difference  equations  in  the 
computational  plane-(ct.P)  for  the  corrector  step  : 
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where 


.  n+1  n+l  n 
Aq  =  u  - u 


A  n+)  n+1  n 

Aqy  =  vy  -  vy 

a,x,  P,t,  p,z  and  P,zz  denote  the  derivatives  of  the  corresponding 
mapping  functions.  The  difference  operators  used  :  V  .  v  <se. 
defined  in  Warming  and  Beam  [34]  : 


8  =  (  )i+l/2  -  (  )i-l/2  5+ =(■),+  ,- (-)i  5' =  (  ),- (  ),., 

52  =  (-)i+i  -  2(  )i  +  (  ),., 


where  x,,  x2  and  D,.  Ds  are  constants,  e  is  a  small  number  which 
controls  the  stretching  in  the  streamwise  coordinate  x.  The  mapping 
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H-  =  1/2  [  (')i+i/2  +  i  h  \r  1 


Equations  (8)  and  (9)  ate  then  factored  into  two  sweeps  by  the 
approximate  factorization  procedure  of  Warming  and  Beam  1 34 1 
Each  sweep  involves  solving  a  system  ot  atgebraic  equations  wuth  a 
scalar  tri-diagonal  coefficient  matrix  : 

ALPlz  g2  _  (2p5)[j)n*1/“  Aq11*1  =  RHS  of  (8) 
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Having  solved  u  and  vy,  the  transverse  velocity  w  can  be  updated 
from  integrating  the  continuity  equation  (1)  by  trapezoidal  rule. 

Similar  difference  schemes  are  applied  to  the  equations  in  the 
transformed  coordinates  systems  for  the  nose  and  tail  regions.  Their 
final  forms  are  omitted  here  for  brevity. 

The  accuracy  of  the  predictor-corrector  method  used  has  been 
studied  by  Douglas  133]  on  1-D  model  equations.  The  results  show 
that  for  a  transient  problem,  the  overall  accuracy  in  time  is  of  second 
order.  Linearized  von  Neumann  stability  analysis  of  the  corrector  step 
predicts  the  scheme  to  be  unconditionally  stable  if  the  one-side 
difference  conforms  with  the  flow  direction.  In  the  following 
numerical  examples,  no  instability  was  encountered  in  the 
calculations.  Furthermore,  in  order  to  preserve  second  order  accuracy 
in  the  spatial  direction,  the  convection  terms  were  discretized  by  a 
three-point  upwind  difference  formula.  This  will  result  in  a  penta- 
diagonal  matrix  in  the  ct-sweep  rather  than  tri-diagonal  during  the 
factorization  procedure.  Both  matrices  w-ere  solved  using  well 
vectorized  penta-  and  tri-diagonal  solvers. 


RESULTS  AND  DISCUSSION 


Calculations  were  done  for  the  symmetry  plane  of  a  prolate 
spheroid  with  axes  ratio  1/4  and  various  angles  of  attack  which  typify 
the  axisymmecry  (0°),  low  incidence  (6°)  and  high  incidence  (50°) 
categories.  To  assure  the  convergence  of  the  solution,  a  numerical  test 
on  the  effects  of  grid  refinement  was  performed  first  and  the  results 
are  included  in  the  appendix.  At  each  different  incidence,  the  skin 
friction  and  displacement  thickness  are  presented  and  discussed  in  the 
following  subsections.  The  focus  will  be  on  the  initiation  of 
separation  evidenced  by  the  emergence  of  a  Lagrangian  singularity  in 
the  flow.  Because  of  the  similarity  between  the  flow  in  the  symmetry 
plane  and  the  two-dimensional  cylinder  case,  comparison  with  the 
results  obtained  from  the  two-dimensional  elliptic  cylinder  with  the 
same  axes  ratio  is  made,  and  certain  analogy  between  the  two  flows  is 
shown. 

All  the  numerical  computations  were  performed  on  the  Cornell 
National  Supercomputer  Facility.  For  a  typical  30,000  grid  points 
(including  all  four  regions)  used  in  each  of  the  present  examples,  the 
required  CPU  time  in  advancing  one  time  step  was  around  2  seconds 
on  IBM  3090-600E.  About  45%  of  the  job  was  executed  under  the 
vector  mode,  and  the  measured  vector  speed-up  was  3.  More  effort  is 
being  made  to  improve  the  percentage  of  vectorization  of  the  program. 


(TJiKidgrtce 

The  zero  degree  incidence  case  was  not  included  ir.  Xu  and 
Wang  s  calculation  [151.  It  is  given  here  as  a  bench  test  of  the  present 
numerical  scheme.  Calculated  wall  shear  and  displacement  thickness 

are  shown  in  Figs. 3  to  8.  In  the  nose  region,  the  skin  friction  i* 
quickly  approaches  a  steady  value  (Fig. 3)  while  in  the  tail  region. 
Fig. 4,  a  reversed-flow  region  is  seen  to  form  shortly  after  the 
impulsive  start  and  spread  subsequently  toward  the  leading  edge. 
Also  notice  that  due  to  the  small  axes  ratio,  the  reversed  flow  is 
confined  in  a  narrow  region  near  the  trailing  edge,  and  the  skin  friction 
does  not  vary  much  along  most  part  of  the  body  (Fig. 5).  Figure  6 
plots  the  distribution  of  crosswise  skin  friction  xyiy.  Owing  to  the 
axisymmetry  property  of  the  flow  at  0°  incidence,  the  value  of  xy,y  is 
zero  in  this  particular  case.  However,  as  pointed  out  in  the  previous 
section,  quantities  on  the  left  and  right  interfaces  of  the  main  region 


are  obtained  from  interpolating  the  flow  variables  in  the  nose  and  tail 
regions  respectively.  Specifically,  the  crosswise  skin  friction  in  the 
main  region  is  related  to  that  in  the  nose  through  the  transformation 

Ty  y  _  -T-  COS  y  +  Ty  y  5  <.OS~  }’. 

In  general,  terms  on  the  right  hand  side  of  the  above  transformation  do 
not  sum  up  to  zero  due  to  discretization  and  hence  introducing  error  on 
the  interface.  It  can  be  seen  from  Fig. 6  that  this  error  on  the  left 
interface  is  quickly  damped  out  after  a  few  nodal  poinrs.  while  the 
error  on  the  downstream  interface  does  not  propagate  back  into  the 
domain  because  of  the  upwind  differencing.  The  accuracy  of  the 
solution  at  the  interior  points  is  found  to  be  essentially  unaffected. 

The  displacement  thickness  AK  is  given  in  Fig. 7  and  Fig. 8.  In  the 
main  region  (Fig.7),  the  flow  changes  slowly  over  a  large  pan  of  the 
body  except  very  close  to  the  trailing  edge  where  the  boundary-layer 
increases  drastically.  Figure  8  shows  the  details  just  ahead  of  the 
trailing  edge.  A  "iKike"-like  structure  similar  to  that  discovered  in  the 
Lagrangian  computations  of  two-dimensional  problems  ([17],  [18], 

[351)  is  observed  near  x=-0.08.  In  order  to  be  convinced  that  it  is  not 
created  by  numerical  instability,  the  calculation  was  carried  out  again 
in  several  refined  meshes  and  reduced  time  steps.  Results  show  only 
a  "sharper''  spike  and  the  overall  picture  does  not  change.  Formerly, 
this  "spike"-like  behavior  has  only  been  found  in  Lagrangian 
computations.  Most  of  the  previous  two-dimensional  boundary-layer 
results  calculated  by  the  Eulerian  formulation  encountered  difficulty 
near  "separation"  and  consequently  did  not  give  a  conclusive 
displacement  thickness  variation  at  this  point  (e.g.  [36-39]).  Failure 
to  obtain  a  converged  and  stable  solution  in  these  calculations  upon 
approaching  separation  may  be  attributed  to  the  aforementioned  severe 
limitation  on  the  step  size  needed  in  the  conventional  space-marching 
technique.  The  present  time-marching  Eulerian  scheme  proves  to  be 
not  only  stable  but  also  capable  of  confirming  the  critical  features  near 
separation  w-hich  were  previously  obtained  by  the  Lagrangian 
computation. 

The  behavior  of  the  displacement  thickness  reflects  the  behavior  of 
the  velocity  profile  u.  Shown  in  Fig.9  is  the  velocity  profile  at  several 
x  locations.  It  is  found  that  starting  from  t=0.19,  there  is  an  unusual 

"clustering"  of  the  velocity  profiles  around  x  =-0.045.  The  same 
phenomenon  was  reported  in  Cebeci  [36].  The  clustering  of  the 

velocity  profiles  results  in  a  "kink"  in  the  Aj  distribution  at  t=0.20. 
Fig. 8.  Beyond  t=0.20,  there  is  a  sudden  thickening  in  some  of  the 
velocity  profiles  (Figs.9b,  c)  which  in  turn  leads  to  a  fast  increase  in 
the  displacement  thickness  at  such  locations.  The  peculiar  velocity 

profiles  are  only  found  to  occur  in  a  very  narrow  region  of  x.  Further 
into  the  reversed  flow  regime,  the  velocity  profiles  look  normal  again. 
The  formation  of  the  sharp  "spike"  and  the  "valley"  immediately 
following  in  the  displacement  thicknesses  for  t>0.20,  Fig.8,  are  thus 
explained. 

Accompanying  this  thickening  of  the  displacement  thickness  is  the 

sudden  growth  in  the  transverse  velocity  w.  From  the  Lagrangian 
view  point  ([18],  [40]),  the  thickening  happens  because  a  certain  fluid 
"packet"  gets  squeezed  in  the  streamwise  direction  and  must  expand 
laterally  in  order  to  preserve  the  volume.  Since  in  the  axisymmetry 
case  every  meridian  plane  that  cuts  through  the  prolate  spheroid  is  a 
symmetry  plane,  the  critical  fluid  packets  form  a  front  that  has  no 
choice  but  to  explode  in  the  direction  normal  to  the  wall.  The 
projection  of  this  front  onto  the  wall  encircles  the  spheroid,  and  is  the 
logical  generalization  of  the  separation  point  for  two-dimensional 
flows,  it  unequivocally  defines  the  "separation  line"  for  unsteady 
three-dimensional  separation. 

The  singular  behavior  of  the  boundary-layer  at  separation  is  also 

seen  from  the  continuously  increasing  3u/3x  near  separation. 
According  to  the  two-dimensional  Lagrangian  computation  of  van 

Dommeien  [18],  3u/3x  should  blow  up  at  separation.  Moreover,  the 
velocity  profile  near  separation  should  exhibit  a  very  "fiat"  front 
corresponding  to  a  large  inviscid  region  that  splits  the  boundary-layer 
into  two  vorticity  layers.  The  present  Eulerian  calculation  did  find  a 
thickening  of  some  velocity  profiles  as  was  shown  in  Fig.9,  but  none 
of  the  profiles  are  very  "flat".  Now,  given  a  grid  size,  a  simple 
estimation  based  on  the  Taylor  series  expansion  would  suggest  an 
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upper  bound  of  the  value  3u/Dx  that  can  be  determined  by 

differencing.  Specifically,  when  the  term  (Au/Ax )( Ax)  becomes  order 
1,  the  finite  difference  formula  ceases  to  be  meaningful.  In  other 
words,  if  the  mesh  Ax  =0.007,  which  was  the  grid  size  near  separation 
in  this  case,  the  maximum  Au/Ax  we  can  get  is  roughly  150.  Indeed 

in  the  calculauon,  it  is  found  that  the  maximum  value  of  Au/Ax  hardly 
exceeded  this  order.  As  such  a  maximum  is  reached,  the  singularity  is 
smeared  out  by  the  numerical  diifusion,  and  the  legitimacy  of 
continuing  the  calculation  becomes  questionable.  Higher  resolution 
thus  is  needed  to  obtain  more  accurate  velocity  profiles  near 
separation.  Further  discussion  about  this  will  be  given  in  the 
subsequent  sections. 

6°  Incidence 

This  can  be  considered  as  an  axisymmetry  flow  perturbed  by  a 
small  angle  of  attack.  It  is  expected  that  most  of  the  features  will  be 
similar  to  those  in  the  0°  incidence.  Figure  10  shows  the  lee-side  skin 
friction,  also  superimposed  is  the  result  by  Xu  and  NVang  [15],  The 
two  results  do  agree  in  the  general  trend  but  not  in  the  values.  The 
skin  friction  in  the  tail  region  is  plotted  in  Fig.l  1  where  details  of  a 
small  region  of  reversed  flow  near  the  trailing  edge  are  revealed. 
(This  was  not  available  in  [15].)  Again,  no  anomalous  phenomenon 
is  found  in  the  wall  shear,  and  the  whole  picture  looks  almost  the 
same  with  the  one  in  the  previous  case  (Fig.4)  except  now  the  rear 
stagnation  point  has  shifted  upstream  a  little  bit  from  the  rear  end.  In 
the  present  numerical  scheme,  the  velocity  profile  at  the  stagnation 
point  is  not  needed  to  construct  the  flow  as  normally  done  in  a  space¬ 
marching  technique.  The  stagnation  flow  evolves  by  itself.  From 
Fia.l  1,  it  is  seen  that  both  the  stagnation  point  and  the  asymmetry  of 
the  flow  due  to  the  small  angle  of  attack  are  well  captured  in  the 
calculation. 

Figure  12  exhibits  a  non-zero  distribution  of  crosswise  skin 
friction  ty  y  at  lee-side.  At  small  times,  cy,y  remains  negative  which 
indicates  that  the  circumferential  velocity  v  near  the  symmetry  plane  is 
pointing  toward  the  plane  of  symmetry.  As  time  increases,  a  region  of 

positive  ty,y  begins  to  develop  near  the  trailing  edge.  In  other  words, 
a  small  crosswise  reversed  flow  is  formed  at  the  rear  end.  This 
reversed  flow  region  will  then  extend  gradually  toward  the  leading 
edge  direction  with  increasing  time.  The  present  results  are  consistent 
with  those  described  in  [15]. 

Both  the  lee-side  and  wind-side  displacement  thicknesses  in  the 
tail  region  are  shown  in  Figs.  13  and  14.  (Details  of  the  displacement 
thickness  near  the  rear  end  were  left  out  in  Xu  and  Wang  [15], 
therefore  no  comparison  is  made.)  By  examining  Fig.  13,  it  is 
observed  that  a  "spike"  appears  first  in  the  wind-side  displacement 
thickness  at  t=0.25.  If  we  accept  the  assumption  that  the  initial 
separation  does  not  substantially  alter  the  flow  over  the  lee-side  and 
simply  proceed  with  the  calculation,  a  similar  "spike"-like  structure  is 
also  found  to  occur  on  the  lee-side  at  a  somewhat  later  time  (Fig.  14). 
Without  the  benefit  of  off-symmetry-plane  solution,  it  may  be 
conjectured  here  that  the  separation  is  likely  to  occur  first  as  a  point  on 
the  wind-side  symmetry  plane  and  then  expand  with  increasing  time  in 
both  circumferential  directions,  and  finally,  meet  on  the  lee-side 
symmetry  plane  to  form  a  closed  curve.  Turning  to  the  Lagrangian 
interpretation,  because  of  the  presence  of  a  small  incidence,  fluid 
"packets"  on  each  different  meridian  plane  will  now  deform 
differently.  The  time  at  which  these  "packets"  might  get  squeezed 
into  zero  thickness  would  also  differ.  Besides,  the  three-dimensional 
geometry  offers  the  fluid  "packets"  a  chance  to  expand  in  the 
crosswise  direction.  The  details  of  how  the  separation,  after  its 
initiation  from  the  wind-side  symmetry  plane,  propagates  along  the 
circumferential  direction  is  a  question  that  can  only  be  answered  after  a 
fully  three-dim-nsional  calculation  over  the  entire  prolate  spheroid  is 
made. 

KT  Incidence 

After  the  impulsive  start-up,  the  lee-side  skin  friction  decreases 
from  a  maximum  value  at  the  leading  edge  to  zero  at  the  rear 
stagnation  point  (Fig. 15).  However,  beginning  from  t=0.05,  a  dip 


near  the  front  end  is  found,  and  the  value  of  tx  at  the  dip  drops 
quickly  thereafter.  According  to  the  large  time  calculation  of  Xu  and 
Wang  [15],  this  value  will  eventually  become  negative  as  a  small 
recirculating  bubble  is  established.  A  similar  situation  was  found  in  a 
high-incidence  elliptic  cylinder  [35]  where  a  reversed  flow  appears  at 
early  time  near  the  leading  edge.  The  skin  friction  in  the  tail  region  is 
presented  in  Fig.  16.  In  this  picture,  the  rear  stagnation  point  is 

located  at  x=-0.23.  By  the  definition  of  the  coordinate  transformation, 
positive  x  points  to  wind-side  direction.  Therefore,  positive  wall 

shear  around  x=0  in  this  plot  actually  corresponds  to  a  recirculating 
flow  near  the  rear  end.  The  recirculation  region  first  occurs  just 
shortly  above  the  trailing  edge  and  then  expand  in  both  directions. 

The  zero  point  moving  toward  the  wind-side  direction  appears  to 

stop  at  x=0.0l,  and  the  recirculation  region  does  not  penetrate  further 
upstream.  In  fact,  the  skin  friction  at  the  wind-side  seems  to  have 
aireadv  reached  a  steady  value  at  t=0.12.  The  same  trend  was  also 
found  in  the  previous  two-dimensional  calculation  [35]. 

An  interesting  behavior  in  the  lee-side  crosswise  skin  friction  is 
demonstrated  in  Fig.  17.  The  reversal  of  circumferential  velocity  v 
(positive  Ty  y)  first  occurs  at  the  rear  end  and  the  zero  ty  y  point 
quickly  moves  upstream.  At  t=0.08,  another  region  of  reversed  v 
appears  at  the  front  and  grows  in  the  downstream  direction.  These 
two  regions  coalesce  at  approximately  t=0.12  which  means  the  flow 
near  the  symmetry  plane  is  moving  away  from  the  plane  over  the 
entire  lee-side  except  for  a  small  part  close  to  the  leading  edge.  It  is 
believed  that  the  reversal  of  the  crosswise  velocity  v  will  have  an 
important  influence  on  the  off-symmetry-plane  separation. 

At  such  a  high  angle  of  attack,  the  lee-side  displacement  thickness 
at  small  time  still  remains  essentially  constant  along  most  part  of  the 
body  (Fig.  18).  Yet  a  "hump"  resulting  from  the  "dip"  in  the  skin 
friction  (Fig.  15)  is  evolved  at  t=0.05.  Near  the  rear  end,  the 
displacement  thickness  shows  a  sharp  increase  and  forms  a  "spike" 
almost  right  at  the  trailing  edge  (Fig.  19).  This  fact  suggests  a  possible 
difficulty  of  the  method  used  by  Xu  and  Wang  [15]  in  obtaining 
accurate  results  near  the  trailing  edge,  because  they  did  not  remove  the 
geometrical  singularity  there.  The  reason  that  the  "spike"  at  t=0. 12 
appears  to  be  sharper  than  those  presented  in  the  previous  cases  of  this 
paper  is  due  to  a  somewhat  higher  resolution  in  the  present  case.  The 

grid  size  Ax  near  the  trailing  edge  here  is  approximately  0.00 1. 
According  to  the  argument  stated  at  the  end  of  the  First  subsection,  the 

maximum  Au/Ax  that  can  be  obtained  here  is  around  1000,  which  is 
an  order  larger  than  that  in  the  zero  incidence  case.  Associated  with 
this  result  is  a  "flatter"  velocity  profile  near  separation  (Figs.20a,  b). 
Although  the  relevant  profiles  are  not  too  smooth,  the  trend  is  clear. 
Shown  in  Fig.21  is  one  of  the  corresponding  vorticity  profiles.  The 
splitting  of  the  boundary-layer  into  two  vorticity  layers  is  apparent. 

The  calculation  was  terminated  at  t=0. 12  after  the  presence  of  a 
singularity  on  the  wind-side.  Whether  or  not  another  singularity  may 
appear  at  the  lee-side  near  the  leading  edge  can  not  be  answered  at  this 
moment.  However,  suggested  by  the  large  time  calculation  of  Xu  and 
Wang  [15]  and  the  former  results  on  the  elliptic  cylinder  [35],  it  is 
likely  that  another  singularity  would  develop  at  later  time  near  the  front 
end,  and  afterwards  the  separation  would  probably  be  of  the  "closed" 
type  at  this  high  incidence. 

Comparison  with  Two-Dimensional  Results 

The  results  of  an  unsteady  boundary-layer  over  an  impulsively 
started  elliptic  cylinder  with  axes  ratio  1/4,  calculated  by  the 
Lagrangian  method  described  in  [35],  are  presented  here  for 
comparison.  Shown  in  Figs. 22a,  b,  c  are  the  displacement 
thicknesses  at  0°,  6°  and  50°  incidence.  In  these  plots,  the  time  t  has 
been  convened  to  conform  with  the  present  time  scale,  %  denotes  the 
streamwise  boundary-layer  coordinate  measured  along  the  body  from 
the  leading  edge  in  the  clockwise  direction.  Resemblance  beov. .-  ■>  T< 
two  cases  is  clear  (cf.  Figs. 8,  14,  19)  except  that  in  the  Lagrantp.m 
result,  the  "spike"  is  steeper  and  more  pronounced.  This  is  i  '  6 
from  the  high  resolution  of  the  Lagrangian  grid  near  separation  <  'q). 
Besides,  the  displacement  thickness  in  the  Lagrangian  calculation  is 
obtained  from  integrating  the  continuity  equation  which  turns  singular 
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as  the  separation  is  approached. 

It  is  noticed  that  tor  all  three  incidences  calculated  here,  separation 
occurs  earlier  in  time  in  the  two-dimensional  case.  This  can  be 
attributed  to  the  curvature  effect  as  well  as  the  non-zero  crosswise 
velocity  vy  in  the  present  three-dimensional  calculation. 

Figure  23  plots  the  maximum  3u/(hdx)  in  the  present  symmetry 
plane  calculation  against  (rs-t),  ts  being  the  separation  time  when  the 

solution  turns  singular.  At  ts,  the  maximum  3u/(h3x)  goes  to  infinity 
theoretically.  We  obtain  an  approximate  ts  by  extrapolating  the  zero 

point  of  l/(<3u/(hdx))  from  the  available  data.  In  the  50°  incidence 
case,  least-square  fittings  with  models  from  quadratic  to  5th  order 
polynomial  give  almost  the  same  value  of  ts  around  0.123.  The  result 
for  the  6°  incidence  is  less  certain  because  of  the  scattering  and 

relatively  low  value  of  maximum  3u/(h3x),  only  a  quadratic 
polynomial  was  used  to  obtain  ts  Also  shown  in  Fig. 28  are  the 
straightlines  of  slope  -7/4  predicted  by  van  Dommelen  (181  from  the 
asymptotic  structure  of  the  two-dimensional  singularity.  The  dotted 

line  indicates  the  limiting  value  of  3u/(hdx)  beyond  which  the 

calculated  result  is  unreliable.  The  blow-up  of  maximum  du/(hdx) 
appears  to  follow  the  7/4  line.  This  suggests  that  the  structure  of  the 
separation  on  the  symmetry  plane  is  similar  to  that  for  the  two- 
dimensional  case. 

SUMMARY 

The  subject  of  unsteady  three-dimensional  boundary-layer 
separation  is  of  grow  ing  importance.  Recent  understanding  about  the 
unsteady  separation  mechanism  reveals  that  the  usual  waif  shear  and 
limiting  streamlines  analysis  are  not  strictly  applicable  to  the  separation 
in  an  unsteady  flow.  An  accurate  and  efficient  computation  scheme 
for  a  general  three-dimensional  boundary-layer  flow  is  needed.  The 
present  paper  is  a  progress  report  of  a  broader  research  program  for 
this  difficult  problem. 

We  propose  a  non-iterative,  time-advancing  numerical  scheme 
along  with  a  multiple-zone  strategy  to  solve  the  unsteady  three- 
dimensional  boundary- layer  equations.  As  a  preliminary  study,  the 
problem  of  Xu  and  Wang  [151,  i  e.  the  symmetry-plane  solution  of  an 
impulsively  started  prolate  spheroid  with  axes  ratio  1/4  has  been 
redone.  Three  different  incidences  0°,  6°  and  50°  are  calculated. 
Results  from  the  examples  indicate  that  the  algorithm  produce  accurate 
results  over  the  enure  symmetrical  plane,  including  certain  Lagrangian 
features  near  separation  found  in  the  two-dimensional  case.  For  the 
0°  incidence,  projection  of  the  separation  locus  off  the  wall  on  the 
body  necessarily  forms  a  closed  curve  encircling  the  spheroid  owing 
to  the  axisymmetry.  For  non-zero  incidence,  separation  could 
conceivably  start  as  a  point  on  the  wind-side  symmetry  plane  and  then 
propagate  toward  the  lee-side  in  the  crossflow  direcnons.  Whether  or 
not  the  curve  will  become  closed  on  reaching  the  lee-side  symmetry 
piune  may  depend  on  the  angle  of  attack.  Based  on  the  behavior  of  the 

maximum  du/(hdx),  the  structure  of  the  three-dimensional  singularity 
on  the  symmetry  plane  seems  to  be  similar  to  that  for  the  two- 
dimensional  case. 

By  using  the  time-advancing  scheme  and  a  flow  dependent  one¬ 
sided  spatial  difference,  the  CFL  condition  is  automatically  satisfied. 
No  instability  was  experienced  in  the  calculations.  The  accuracy  of 
the  solution  near  separation  is  found  to  be  limited  by  the  requirement 
of  increasing  spatial  resolution,  not  by  the  stability  consideration  as 
was  inherent  in  a  space-marching  technique.  Besides,  the  present 
time-advancing  scheme  allows  more  general  flow  motions  to  be 
calculated  while  the  conventional  space-marching  scheme  has  the 
difficulty  of  defining  a  starting  profile  to  initiate  the  spatial  integration. 

The  primary  concern  of  adopting  a  multi-zone  approach  here  is  to 
remove  the  geometrical  singularity  presented  on  the  prolate  spheroid. 
Since  the  calculation  task  in  each  region  is  essentially  independent, 
possible  implementation  of  a  parallel  processing  on  these  multiple 
regions  is  therefore  worth  studying. 
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APPENDIX 

In  this  appendix,  accuracy  of  the  present  numerical  calculations  is 
estimated  by  studying  the  effect  of  grid  refinement  on  the  value  of  -Am 
friction  at  the  front  stagnation  point.  On  the  stagnation  line,  the 
velocity  component  u  is  identically  zero.  The  derivative  of  ttv  skip, 
friction  at  stagnation  is  obtained  by 

fa  =  (I*  *  /  h)  /  (Ue  *  /h  r'2. 

The  scaling  factor  (Ue  y  /  h  )j/2  is  introduced  here  for  the  purpose  of 
comparison  with  Howarth's  steady  state  resuii  [41],  At  steady  state. 

tq  is  equivalent  to  f(0)  defined  in  [4 1  ].  For  the  axisvmmetnc.il  case, 
f '(01=  1.3 12. 

Table  1  lists  x0  at  t=0.2  calculated  at  three  different  mesh  sizes 
with  time  step  At  equals  0.02,  0.01  and  0.005  respectively.  All  three 
values  appear  to  approach  the  steady  state  result  1.312  given  by 
Howarth.  The  last  row  of  the  table  is  the  Richardson  extrapolation  of 
the  above  three  meshes.  Using  this  value  as  a  reference  "exact" 

solution,  errors  in  t0  obtained  from  each  gnd  size  are  listed  in  the 
secon^  column.  Convergence  of  the  result  is  apparent. 

Table  1.  Effects  of  grid  refinement  on  the  value  of  skin  friction  at 
stagnation  point. 


gnd  size  V,  error 

30  x  31  1.310660  0.001269 

60  x  61  1.311750  0.000179 

120  x121  1.311893  0.000036 

RE  1.311929 


Fig.l  Multiple  zones. 
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I .  Introduction 


Since  the  noniterative  factorization  method  for  the  compressible 
Navier -Stokes  equation  was  developed  by  Beam  &  Warmming(l),  several 
similar  noniterative  methods  have  been  investigated  to  solve  the 
compressible  and  the  incompressible  flows.  In  Beam  &  Warmming(l) 1 s 
technique,  the  governing  equation  can  be  noniteratively  solved  without 
loss  of  accuracy.  This  scheme  is  a  kind  of  three  point  backward  implicit 
ADI  method  developed  by  McDonald  and  Briley(2). 

But,  even  though  the  Navier-Stokes  can  be  solved  without 
iterations,  a  large  amount  of  computer  time  and  storage  is  needed 
to  solve  any  simple  flow.  Though  only  the  steady  solution  is  needed, 
the  unsteady  equation  is  calculated  in  time  until  a  steady  state  solution 
is  achieved.  Moreover,  it  is  possible  to  obtain  the  accurate  solutions 
of  many  viscous  flows  with  several  reduced  sets  of  equations,  such 
as  boundary  layer  equations,  thin-layer,  or  parabolized  Navier-Stokes 
equations.  So,  many  researchers  have  developed  the  noniterative  methods 
to  solve  the  reduced  sets  of  equations. 

These  noniterative  methods,  which  have  been  developed  on  the 
base  of  Beam  &  Warmming’s  delta  scheme,  can  be  divided  into  two 
catagories,  i.  e.  one  for  compressible  flows,  and  another  for 
incompressible  flows.  Steger(3)  used  the  technique  to  solve  the  unsteady 
thin-layer  Navier-Stokes  equation.  The  thin-layer  Navier-Stokes  equation 
is  considerably  less  complicated  than  the  complete  Navier-Stokes 
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equation  because  of  the  approximation  of  dropping  some  viscous  terms  in 
the  Navier-Stokes  equation.  But,  the  procedure  of  calculation  still 
requires  a  substantial  amount  of  computer  effort  to  solve  three 
dimensional  flows.  So,  Schiff  and  Steger(4)  developed  the  noniterative 
scheme  to  solve  the  parabolized  Navier-Stokes  equation  to  predict  three 
dimensional,  steady,  supersonic  viscous  flowfield.  This  method  requires 
that  the  inviscid  outer  flow  must  be  supersonic  and  the  primary 
streamwise  velocity  component  must  be  positive  everywhere.  But  the  cross 
flow  separation  is  permitted.  An  additional  constraint  is  that  the 
streamwise  pressure  gradient  in  the  region  of  subsonic  must  be  handled 
to  neglect  the  elliptic  characteristics  of  the  governing  equation. 

As  above,  though  many  noniterative  schemes  have  been  developed  and 
popularly  used  to  solve  compressible  flows,  only  a  few  schemes  have  been 
developled  for  incompressible  flows.  Kim  and  Chang(5)  used  the 
noniterative  technique  to  solve  for  unsteady  incompressible  boundary 
layers.  But,  it  can  not  be  applied  to  the  flow  which  includes  the  region 
of  reverse  flow,  exept  for  some  spetial  cases.  Kwak(6)  developed  the 
technique  to  solve  the  steady  incompressible  Navier-Stokes  equation 
using  the  artificial  compressibility  method  of  Chorin(7).  Although  this 
method  used  the  linearized  factorization  technique,  it  is  the  iterative 
scheme  to  obtain  steady  solutions. 

Some  iterative  methods  have  been  developed  to  solve  the  steady 
parabolized  Navier-Stokes  equation.  The  streamwise  derivative  term  of 
pressure  is  treated  explicitly  in  these  methods.  In  Patankar  and 
Spalding(9),  the  pressure  in  the  streamwise  momentum  equation  is 
assumed  to  vary  only  in  the  streamwise  direction.  But,  the  variation  of 
pressure  in  a  cross  plane  is  permitted.  Chilukuri  and  Pletcher(lO) 
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considered  Che  effect  of  the  elliptic  characteristic  by  marching  in  the 
streamwise  direction  in  an  iterative  method.  This  scheme  was  originally 
restricted  to  flows  in  which  the  reverse  flow  in  the  primary  flow 
direction  doesnot  occur.  But,  Madavan  and  Pletcher(ll)  have  demonstrated 
two  dimensional  flows  in  which  the  reverse  flow  in  the  primary  flow 
direction  occurs.  But,  all  of  these  methods  are  iterative  methods  for 
steady  solutions. 

Here,  a  code  is  developed  to  solve  the  unsteady  3-dimensional 
thin-layer  Navier-Stokes  equation  by  using  the  noniterative  finite 
difference  method.  At  first,  the  Navier-Stokes  equation  written  in  the 
general  coordinates  <,)  is  simplified  by  the  thin-layer 

approximation,  which  neglect  the  viscous  terms  including  the  derivatives 
of  the  streamwise  or  the  spanwise  direction.  This  thin- layer  equation  is 
linearized  and  factorized  by  the  noniterative  technique.  The  resulting 
factorized  equation  is  solved  by  a  3-step  block  tri-diagonal  matrix 
elimination.  The  central  difference  is  used  to  discretize  the 
derivatives  of  the  spanwise  and  the  normal  directions,  but  the  central 
or  backward  difference  is  used  for  the  streamwise  direction. 

The  flow  which  has  a  leading  edge  stagnation  point  can  not  be  solved  by 
the  central  difference,  because  the  results  are  unstable  between  the 
initial  point  and  the  leading  edge  stagnation  point.  So,  the  backward 
difference  is  applied  to  the  derivatives  of  velocities  in  the 
momentum  equations  between  the  initial  point  and  the  leading  edge 
stagnation  point.  This  method  helps  to  stabilize  the  scheme,  but  the 
results  are  slowly  disturbed  again  from  the  leading  edge  stagnation 
point . 
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Some  numerical  examples  are  investigated,  such  as  the  transient 
development  of  Couette  flow,  the  flow  over  a  suddenly  started  finite 
flat  plate,  the  flow  over  a  suddenly  started  circular  cylinder  and 
the  transient  flow  ever  a  3-dimensional  wave  wall. 
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